﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using AdaptiveDeNet.RAND;
using AdaptiveDeNet.CSV;
using System.Threading;
using System.Threading.Tasks;
using System.Diagnostics;

namespace AdaptiveDeNet.DifferentialEvolution.MultiThreaded
{
    public partial class Optimizer
    {
        private object locker = new object();

        public CEC14Result JADE(double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double c, double p, Queue<double> FESMultipliers, double fX, int funcNum)
        {
            mre = new ManualResetEvent(false);
            countdown = new CountdownEvent(threadCount);
            outerCountdown = new CountdownEvent(threadCount);

            lockers = new ManualResetEventSlim[threadCount];

            mreslis = new ManualResetEventSlim[threadCount];

            condition = true;

            double uCr = 0.5;
            double uF = 0.5;

            CEC14Result result = new CEC14Result();

            for (int i = 0; i < mreslis.Length; i++)
            {
                mreslis[i] = new ManualResetEventSlim(false);
            }

            RngStream stream = new RngStream();
            long seed = DateTime.Now.Ticks & Int32.MaxValue;
            RngStream.setPackageSeed(new long[] { seed, seed, seed, seed, seed, seed });

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            double[] uBlBDiff = new double[dim];

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            int popEvals = N;

            //mainPopHolder
            int padding = 16;
            int spacing = 16;

            //pointers
            int[] populationPointers = new int[N];
            int[] fitnessPointers = new int[N];
            int[] SArchivePointers = new int[N];

            double[] mainPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] trialsPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] A = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            int sizeA = 0;
            double[] Scr = new double[padding + (1 * N) + (threadCount * spacing)];
            double[] Sf = new double[padding + (1 * N) + (threadCount * spacing)];

            int archiveHelper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;

                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    SArchivePointers[archiveHelper] = padding + i + threadActual * spacing;
                    archiveHelper++;
                }
            }

            int helper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;
                double[] temp = new double[dim];
                double fEvalInner = 0;
                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    //member vector initialization
                    int start = padding + i * (dim + 1) + threadActual * spacing;
                    int end = start + dim;

                    populationPointers[helper] = start;
                    fitnessPointers[helper] = end;

                    for (int j = start, forDiff = 0; j < end; j++, forDiff++)
                    {
                        temp[forDiff] = nRand.NextDouble() * uBlBDiff[forDiff] + lowerBound[forDiff];
                    }

                    fEvalInner = fEval.Feval(temp); //evaluation
                    if (fEvalInner < Fmin) Fmin = fEvalInner; //is fMin ?


                    Array.Copy(temp, 0, mainPopHolder, start, dim); // temp -> mainPopHolder
                    mainPopHolder[end] = fEvalInner; // fEval -> mainPopHolder

                    helper++;
                }
            }

            Fmax = this.FindMax(mainPopHolder, fitnessPointers);
            Fmin = this.FindMin(mainPopHolder, fitnessPointers);

            //
            //
            //

            long[] affinity = new long[threadCount];

            affinity[0] = 0x0001;
            affinity[1] = 0x0004;


            double[] ScrHolder = new double[padding + (N) + (threadCount * spacing)];
            double[] SfHolder = new double[padding + (N) + (threadCount * spacing)];

            for (int k = 0; k < threadCount; k++)
            {
                int iThread = k;

                ThreadPool.QueueUserWorkItem((st) =>
                {
                    lockers[iThread] = new ManualResetEventSlim(true);


                    int start = N * iThread / threadCount;
                    int end = N * (iThread + 1) / threadCount;
                    int n1 = N - 1;
                    int dim1 = dim - 1;
                    RngStream localStream = new RngStream();

                    double Cr = 0.0;
                    double F = 0.0;

                    uint unmanagedId = GetCurrentThreadId();

                    ProcessThread myThread = (from ProcessThread entry in Process.GetCurrentProcess().Threads
                                              where entry.Id == unmanagedId
                                              select entry).First();

                    myThread.ProcessorAffinity = (IntPtr)affinity[iThread];

                    CEC14Lib.CEC14func localFunc = new CEC14Lib.CEC14func(30, funcNum);

                    double[] temp = new double[dim];
                    double[] tempA = new double[dim];
                    double innerFEval = 0;
                    int i = start;

                    countdown.Signal();

                    do
                    {
                        lockers[iThread].Wait();
                        lockers[iThread].Reset();

                        for (i = start; i < end; i++)
                        {
                            //individual
                            Cr = MathNet.Numerics.Distributions.Normal.Sample(localStream, uCr, 0.1);
                            F = MathNet.Numerics.Distributions.Cauchy.Sample(localStream, uF, 0.1);

                            if (Cr < 0) Cr = 0.0;
                            if (Cr > 1) Cr = 1.0;

                            while (F <= 0) F = MathNet.Numerics.Distributions.Cauchy.Sample(localStream, uF, 0.1);
                            if (F >= 1) F = 1;

                            //get best Pmember
                            int bestP = GetBestPMemberIndex(mainPopHolder, fitnessPointers, p, N, localStream);
                            int r1 = localStream.Next(N);

                            int PUnionAlength = N + sizeA;

                            int r2 = localStream.Next(PUnionAlength);

                            while (i == r1) r1 = localStream.Next(N);
                            while (i == r1 || r1 == r2) r2 = localStream.Next(PUnionAlength);

                            int pom = r2;

                            bool lieInAset = false;

                            if (r2 > (N - 1))
                            {
                                lieInAset = true;
                                r2 = r2 - (N);
                            }

                            if (lieInAset)
                            {
                                lock (locker)
                                {
                                    for (int l = 0; l < dim; l++)
                                    {
                                        tempA[l] = A[populationPointers[r2] + l];
                                    }
                                }
                            }

                            int jrand = localStream.randInt(0, dim1);

                            int innerStart = padding + i * (dim + 1) + iThread * spacing;
                            int innerEnd = innerStart + dim;

                            int popPtr = i;

                            for (int j = innerStart, jPtr = 0; j < innerEnd; j++, jPtr++)
                            {

                                //krizeni
                                if (localStream.randU01() <= Cr || jPtr == jrand)
                                {
                                    double r2val;

                                    if (lieInAset)
                                    {
                                        r2val = tempA[jPtr];
                                    }
                                    else
                                    {
                                        r2val = mainPopHolder[populationPointers[r2] + jPtr];
                                    }

                                    // x = Rcurrent + F * (Rpbest - Rcurrent) + F * (r1 - r2)
                                    temp[jPtr] = mainPopHolder[populationPointers[i] + jPtr]
                                        + F * (mainPopHolder[populationPointers[bestP] + jPtr]
                                        - mainPopHolder[populationPointers[i] + jPtr])
                                + F * (mainPopHolder[populationPointers[r1] + jPtr] - r2val);

                                    if (temp[jPtr] < lowerBound[jPtr] || temp[jPtr] > upperBound[jPtr])
                                    {
                                        temp[jPtr] = Common.Common.Zrdcad(temp[jPtr],
                                            lowerBound[jPtr], upperBound[jPtr]);
                                    }
                                }
                                else
                                {
                                    temp[jPtr] = mainPopHolder[j];
                                }

                            }

                            innerFEval = localFunc.FEval(temp, 30, funcNum);

                            if (innerFEval <= mainPopHolder[innerEnd])
                            {
                                trialsPopHolder[innerEnd] = innerFEval;
                                Array.Copy(temp, 0, trialsPopHolder, innerStart, dim);

                                // old -> A
                                if (sizeA >= N)
                                {
                                    int tempIndex = localStream.Next(N);
                                    lock (locker)
                                    {
                                        Array.Copy(mainPopHolder, innerStart, A, populationPointers[tempIndex], dim);
                                    }

                                }
                                else
                                {
                                    if (sizeA < N)
                                    {
                                        Array.Copy(mainPopHolder, innerStart, A, populationPointers[sizeA], dim);
                                        Interlocked.Increment(ref sizeA);
                                    }
                                }

                                Scr[SArchivePointers[i]] = Cr;
                                Sf[SArchivePointers[i]] = F;

                            }//succesful trial
                            else
                            {
                                trialsPopHolder[innerEnd] = mainPopHolder[innerEnd];
                                Array.Copy(mainPopHolder, innerStart, trialsPopHolder, innerStart, dim);
                            }

                        }


                        mreslis[iThread].Set();
                        //countdown.Signal();

                    } while (condition);

                    outerCountdown.Signal();

                });
            }

            countdown.Wait();

            while (popEvals < maxEvals)
            {
                for (int i = 0; i < threadCount; i++)
                {
                    lockers[i].Set();
                }

                for (int i = 0; i < mreslis.Length; i++)
                {

                    mreslis[i].Wait();
                    mreslis[i].Reset();
                }

                double[] pom;
                pom = mainPopHolder;
                mainPopHolder = trialsPopHolder;
                trialsPopHolder = pom;

                Fmax = this.FindMax(mainPopHolder, fitnessPointers);
                Fmin = this.FindMin(mainPopHolder, fitnessPointers);

                //uCr - computation
                uCr = (1 - c) * uCr + c * meanA(Scr, popEvals, SArchivePointers, N);
                //uF - computation
                uF = (1 - c) * uF + c * meanL(Sf, SArchivePointers, N);

                popEvals += N;

                if (popEvals == (FESMultipliers.Peek() * maxEvals))
                {
                    result.Records.Add(Fmin - fX);
                    FESMultipliers.Dequeue();
                }

                //countdown.Reset(2);
            }

            condition = false;
            //Monitor.PulseAll(locker);
            //countdown.Wait();

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i].Set();
            }

            outerCountdown.Wait();

            result.FuncEvals = popEvals;


            return result;
        }

        /// <summary>
        /// JADE - Adaptive Differential Evolution with external archive,  tries to find global minimum point of given function
        /// </summary>
        /// <param name="lowerBound">Lower bound vector</param>
        /// <param name="upperBound">Upper bound vector</param>
        /// <param name="maxEvals">maximum evaluations of objective function</param>
        /// <param name="dim">dimension of problem</param>
        /// <param name="N">size of population</param>
        /// <param name="myEpps">small nonnegative value, when (Fmax - Fmin) is smaller then myEpps, the termination criteria are satisfied</param>
        /// <param name="c">controls the rate of parameter adaptation, value between 0 and 1</param>
        /// <param name="p">determines the greediness of the mutation strategy (count of selected bestPmembers), value between 0 and 1<</param>
        /// <returns></returns>
        public Result JADE(double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double c, double p)
        {
            mre = new ManualResetEvent(false);
            countdown = new CountdownEvent(threadCount);
            outerCountdown = new CountdownEvent(threadCount);

            lockers = new ManualResetEventSlim[threadCount];

            mreslis = new ManualResetEventSlim[threadCount];

            condition = true;

            double uCr = 0.5;
            double uF = 0.5;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i] = new ManualResetEventSlim(false);
            }


            for (int i = 0; i < mreslis.Length; i++)
            {
                mreslis[i] = new ManualResetEventSlim(false);
            }

            RngStream stream = new RngStream();
            long seed = DateTime.Now.Ticks & Int32.MaxValue;
            RngStream.setPackageSeed(new long[] { seed, seed, seed, seed, seed, seed });

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            double[] uBlBDiff = new double[dim];

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            int popEvals = N;

            //mainPopHolder
            int padding = 16;
            int spacing = 16;

            //pointers
            int[] populationPointers = new int[N];
            int[] fitnessPointers = new int[N];
            int[] SArchivePointers = new int[N];

            double[] mainPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] trialsPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] A = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            int sizeA = 0;
            double[] Scr = new double[padding + (1 * N) + (threadCount * spacing)];
            double[] Sf = new double[padding + (1 * N) + (threadCount * spacing)];

            int archiveHelper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;

                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    SArchivePointers[archiveHelper] = padding + i + threadActual * spacing;
                    archiveHelper++;
                }
            }

            int helper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;
                double[] temp = new double[dim];
                double fEvalInner = 0;
                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    //member vector initialization
                    int start = padding + i * (dim + 1) + threadActual * spacing;
                    int end = start + dim;

                    populationPointers[helper] = start;
                    fitnessPointers[helper] = end;

                    for (int j = start, forDiff = 0; j < end; j++, forDiff++)
                    {
                        temp[forDiff] = nRand.NextDouble() * uBlBDiff[forDiff] + lowerBound[forDiff];
                    }

                    fEvalInner = fEval.Feval(temp); //evaluation
                    if (fEvalInner < Fmin) Fmin = fEvalInner; //is fMin ?


                    Array.Copy(temp, 0, mainPopHolder, start, dim); // temp -> mainPopHolder
                    mainPopHolder[end] = fEvalInner; // fEval -> mainPopHolder

                    helper++;
                }
            }

            Fmax = this.FindMax(mainPopHolder, fitnessPointers);
            Fmin = this.FindMin(mainPopHolder, fitnessPointers);

            double[] ScrHolder = new double[padding + (N) + (threadCount * spacing)];
            double[] SfHolder = new double[padding + (N) + (threadCount * spacing)];

            for (int k = 0; k < threadCount; k++)
            {
                int iThread = k;

                ThreadPool.QueueUserWorkItem((st) =>
                {
                    lockers[iThread] = new ManualResetEventSlim(true);


                    int start = N * iThread / threadCount;
                    int end = N * (iThread + 1) / threadCount;
                    int n1 = N - 1;
                    int dim1 = dim - 1;
                    RngStream localStream = new RngStream();

                    double Cr = 0.0;
                    double F = 0.0;

                    double[] temp = new double[dim];
                    double[] tempA = new double[dim];
                    double innerFEval = 0;
                    int i = start;

                    countdown.Signal();

                    do
                    {
                        lockers[iThread].Wait();
                        lockers[iThread].Reset();

                        for (i = start; i < end; i++)
                        {
                            //individual
                            Cr = MathNet.Numerics.Distributions.Normal.Sample(localStream, uCr, 0.1);
                            F = MathNet.Numerics.Distributions.Cauchy.Sample(localStream, uF, 0.1);

                            if (Cr < 0) Cr = 0.0;
                            if (Cr > 1) Cr = 1.0;

                            while (F <= 0) F = MathNet.Numerics.Distributions.Cauchy.Sample(localStream, uF, 0.1);
                            if (F >= 1) F = 1;

                            //get best Pmember
                            int bestP = GetBestPMemberIndex(mainPopHolder, fitnessPointers, p, N, localStream);
                            int r1 = localStream.Next(N);

                            int PUnionAlength = N + sizeA;

                            int r2 = localStream.Next(PUnionAlength);

                            while (i == r1) r1 = localStream.Next(N);
                            while (i == r1 || r1 == r2) r2 = localStream.Next(PUnionAlength);

                            int pom = r2;

                            bool lieInAset = false;

                            if (r2 > (N - 1))
                            {
                                lieInAset = true;
                                r2 = r2 - (N);
                            }

                            if (lieInAset)
                            {
                                lock (locker)
                                {
                                    for (int l = 0; l < dim; l++)
                                    {
                                        tempA[l] = A[populationPointers[r2] + l];
                                    }
                                }
                            }

                            int jrand = localStream.randInt(0, dim1);

                            int innerStart = padding + i * (dim + 1) + iThread * spacing;
                            int innerEnd = innerStart + dim;

                            int popPtr = i;

                            for (int j = innerStart, jPtr = 0; j < innerEnd; j++, jPtr++)
                            {

                                //krizeni
                                if (localStream.randU01() <= Cr || jPtr == jrand)
                                {
                                    double r2val;

                                    if (lieInAset)
                                    {
                                        r2val = tempA[jPtr];
                                    }
                                    else
                                    {
                                        r2val = mainPopHolder[populationPointers[r2] + jPtr];
                                    }

                                    // x = Rcurrent + F * (Rpbest - Rcurrent) + F * (r1 - r2)
                                    temp[jPtr] = mainPopHolder[populationPointers[i] + jPtr]
                                        + F * (mainPopHolder[populationPointers[bestP] + jPtr]
                                        - mainPopHolder[populationPointers[i] + jPtr])
                                + F * (mainPopHolder[populationPointers[r1] + jPtr] - r2val);

                                    if (temp[jPtr] < lowerBound[jPtr] || temp[jPtr] > upperBound[jPtr])
                                    {
                                        temp[jPtr] = Common.Common.Zrdcad(temp[jPtr],
                                            lowerBound[jPtr], upperBound[jPtr]);
                                    }
                                }
                                else
                                {
                                    temp[jPtr] = mainPopHolder[j];
                                }

                            }

                            innerFEval = fEval.Feval(temp);

                            if (innerFEval <= mainPopHolder[innerEnd])
                            {
                                trialsPopHolder[innerEnd] = innerFEval;
                                Array.Copy(temp, 0, trialsPopHolder, innerStart, dim);

                                // old -> A
                                if (sizeA >= N)
                                {
                                    int tempIndex = localStream.Next(N);
                                    lock (locker)
                                    {
                                        Array.Copy(mainPopHolder, innerStart, A, populationPointers[tempIndex], dim);
                                    }

                                }
                                else
                                {
                                    if (sizeA < N)
                                    {
                                        Array.Copy(mainPopHolder, innerStart, A, populationPointers[sizeA], dim);
                                        Interlocked.Increment(ref sizeA);
                                    }
                                }

                                Scr[SArchivePointers[i]] = Cr;
                                Sf[SArchivePointers[i]] = F;

                            }//succesful trial
                            else
                            {
                                trialsPopHolder[innerEnd] = mainPopHolder[innerEnd];
                                Array.Copy(mainPopHolder, innerStart, trialsPopHolder, innerStart, dim);
                            }

                        }

                        mreslis[iThread].Set();

                    } while (condition);

                    outerCountdown.Signal();

                });
            }

            countdown.Wait();

            while (popEvals < maxEvals)
            {
                for (int i = 0; i < threadCount; i++)
                {
                    lockers[i].Set();
                }

                for (int i = 0; i < mreslis.Length; i++)
                {

                    mreslis[i].Wait();
                    mreslis[i].Reset();
                }

                double[] pom;
                pom = mainPopHolder;
                mainPopHolder = trialsPopHolder;
                trialsPopHolder = pom;

                Fmax = this.FindMax(mainPopHolder, fitnessPointers);
                Fmin = this.FindMin(mainPopHolder, fitnessPointers);

                //uCr - computation
                uCr = (1 - c) * uCr + c * meanA(Scr, popEvals, SArchivePointers, N);
                //uF - computation
                uF = (1 - c) * uF + c * meanL(Sf, SArchivePointers, N);

                popEvals += N;
            }

            condition = false;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i].Set();
            }

            outerCountdown.Wait();

            Result result = new Result();

            result.FuncEvals = popEvals;
            result.Fmin = Fmin;
            result.Point = Common.Common.GetBestMember(mainPopHolder, fitnessPointers, N, dim);


            return result;
        }
        /// <summary>
        /// JADE - Adaptive Differential Evolution with external archive,  tries to find global minimum point of given function
        /// </summary>
        /// <param name="function">objective function</param>
        /// <param name="lowerBound">Lower bound vector</param>
        /// <param name="upperBound">Upper bound vector</param>
        /// <param name="maxEvals">maximum evaluations of objective function</param>
        /// <param name="dim">dimension of problem</param>
        /// <param name="N">size of population</param>
        /// <param name="myEpps">small nonnegative value, when (Fmax - Fmin) is smaller then myEpps, the termination criteria are satisfied</param>
        /// <param name="c">controls the rate of parameter adaptation, value between 0 and 1</param>
        /// <param name="p">determines the greediness of the mutation strategy (count of selected bestPmembers), value between 0 and 1<</param>
        /// <returns></returns>
        public Result JADE(Func<double[], double> function, double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double c, double p)
        {
            mre = new ManualResetEvent(false);
            countdown = new CountdownEvent(threadCount);
            outerCountdown = new CountdownEvent(threadCount);

            lockers = new ManualResetEventSlim[threadCount];

            mreslis = new ManualResetEventSlim[threadCount];

            condition = true;

            double uCr = 0.5;
            double uF = 0.5;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i] = new ManualResetEventSlim(false);
            }


            for (int i = 0; i < mreslis.Length; i++)
            {
                mreslis[i] = new ManualResetEventSlim(false);
            }


            RngStream stream = new RngStream();
            long seed = DateTime.Now.Ticks & Int32.MaxValue;
            RngStream.setPackageSeed(new long[] { seed, seed, seed, seed, seed, seed });

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            double[] uBlBDiff = new double[dim];

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            int popEvals = N;

            //mainPopHolder
            int padding = 16;
            int spacing = 16;

            //pointers
            int[] populationPointers = new int[N];
            int[] fitnessPointers = new int[N];
            int[] SArchivePointers = new int[N];

            double[] mainPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] trialsPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] A = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            int sizeA = 0;
            double[] Scr = new double[padding + (1 * N) + (threadCount * spacing)];
            double[] Sf = new double[padding + (1 * N) + (threadCount * spacing)];

            int archiveHelper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;

                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    SArchivePointers[archiveHelper] = padding + i + threadActual * spacing;
                    archiveHelper++;
                }
            }

            int helper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;
                double[] temp = new double[dim];
                double fEvalInner = 0;
                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    //member vector initialization
                    int start = padding + i * (dim + 1) + threadActual * spacing;
                    int end = start + dim;

                    populationPointers[helper] = start;
                    fitnessPointers[helper] = end;

                    for (int j = start, forDiff = 0; j < end; j++, forDiff++)
                    {
                        temp[forDiff] = nRand.NextDouble() * uBlBDiff[forDiff] + lowerBound[forDiff];
                    }

                    fEvalInner = function.Invoke(temp); //evaluation
                    if (fEvalInner < Fmin) Fmin = fEvalInner; //is fMin ?


                    Array.Copy(temp, 0, mainPopHolder, start, dim); // temp -> mainPopHolder
                    mainPopHolder[end] = fEvalInner; // fEval -> mainPopHolder

                    helper++;
                }
            }

            Fmax = this.FindMax(mainPopHolder, fitnessPointers);
            Fmin = this.FindMin(mainPopHolder, fitnessPointers);

            double[] ScrHolder = new double[padding + (N) + (threadCount * spacing)];
            double[] SfHolder = new double[padding + (N) + (threadCount * spacing)];

            for (int k = 0; k < threadCount; k++)
            {
                int iThread = k;

                ThreadPool.QueueUserWorkItem((st) =>
                {
                    lockers[iThread] = new ManualResetEventSlim(true);


                    int start = N * iThread / threadCount;
                    int end = N * (iThread + 1) / threadCount;
                    int n1 = N - 1;
                    int dim1 = dim - 1;
                    RngStream localStream = new RngStream();

                    double Cr = 0.0;
                    double F = 0.0;

                    double[] temp = new double[dim];
                    double[] tempA = new double[dim];
                    double innerFEval = 0;
                    int i = start;

                    countdown.Signal();

                    do
                    {
                        lockers[iThread].Wait();
                        lockers[iThread].Reset();

                        for (i = start; i < end; i++)
                        {
                            //individual
                            Cr = MathNet.Numerics.Distributions.Normal.Sample(localStream, uCr, 0.1);
                            F = MathNet.Numerics.Distributions.Cauchy.Sample(localStream, uF, 0.1);

                            if (Cr < 0) Cr = 0.0;
                            if (Cr > 1) Cr = 1.0;

                            while (F <= 0) F = MathNet.Numerics.Distributions.Cauchy.Sample(localStream, uF, 0.1);
                            if (F >= 1) F = 1;

                            //get best Pmember
                            int bestP = GetBestPMemberIndex(mainPopHolder, fitnessPointers, p, N, localStream);
                            int r1 = localStream.Next(N);

                            int PUnionAlength = N + sizeA;

                            int r2 = localStream.Next(PUnionAlength);

                            while (i == r1) r1 = localStream.Next(N);
                            while (i == r1 || r1 == r2) r2 = localStream.Next(PUnionAlength);

                            int pom = r2;

                            bool lieInAset = false;

                            if (r2 > (N - 1))
                            {
                                lieInAset = true;
                                r2 = r2 - (N);
                            }

                            if (lieInAset)
                            {
                                lock (locker)
                                {
                                    for (int l = 0; l < dim; l++)
                                    {
                                        tempA[l] = A[populationPointers[r2] + l];
                                    }
                                }
                            }

                            int jrand = localStream.randInt(0, dim1);

                            int innerStart = padding + i * (dim + 1) + iThread * spacing;
                            int innerEnd = innerStart + dim;

                            int popPtr = i;

                            for (int j = innerStart, jPtr = 0; j < innerEnd; j++, jPtr++)
                            {

                                //krizeni
                                if (localStream.randU01() <= Cr || jPtr == jrand)
                                {
                                    double r2val;

                                    if (lieInAset)
                                    {
                                        r2val = tempA[jPtr];
                                    }
                                    else
                                    {
                                        r2val = mainPopHolder[populationPointers[r2] + jPtr];
                                    }

                                    // x = Rcurrent + F * (Rpbest - Rcurrent) + F * (r1 - r2)
                                    temp[jPtr] = mainPopHolder[populationPointers[i] + jPtr]
                                        + F * (mainPopHolder[populationPointers[bestP] + jPtr]
                                        - mainPopHolder[populationPointers[i] + jPtr])
                                + F * (mainPopHolder[populationPointers[r1] + jPtr] - r2val);

                                    if (temp[jPtr] < lowerBound[jPtr] || temp[jPtr] > upperBound[jPtr])
                                    {
                                        temp[jPtr] = Common.Common.Zrdcad(temp[jPtr],
                                            lowerBound[jPtr], upperBound[jPtr]);
                                    }
                                }
                                else
                                {
                                    temp[jPtr] = mainPopHolder[j];
                                }

                            }

                            innerFEval = function.Invoke(temp);

                            if (innerFEval <= mainPopHolder[innerEnd])
                            {
                                trialsPopHolder[innerEnd] = innerFEval;
                                Array.Copy(temp, 0, trialsPopHolder, innerStart, dim);

                                // old -> A
                                if (sizeA >= N)
                                {
                                    int tempIndex = localStream.Next(N);
                                    lock (locker)
                                    {
                                        Array.Copy(mainPopHolder, innerStart, A, populationPointers[tempIndex], dim);
                                    }

                                }
                                else
                                {
                                    if (sizeA < N)
                                    {
                                        Array.Copy(mainPopHolder, innerStart, A, populationPointers[sizeA], dim);
                                        Interlocked.Increment(ref sizeA);
                                    }
                                }

                                Scr[SArchivePointers[i]] = Cr;
                                Sf[SArchivePointers[i]] = F;

                            }//succesful trial
                            else
                            {
                                trialsPopHolder[innerEnd] = mainPopHolder[innerEnd];
                                Array.Copy(mainPopHolder, innerStart, trialsPopHolder, innerStart, dim);
                            }

                        }

                        mreslis[iThread].Set();

                    } while (condition);

                    outerCountdown.Signal();

                });
            }

            countdown.Wait();

            while (popEvals < maxEvals)
            {
                for (int i = 0; i < threadCount; i++)
                {
                    lockers[i].Set();
                }

                for (int i = 0; i < mreslis.Length; i++)
                {

                    mreslis[i].Wait();
                    mreslis[i].Reset();
                }

                double[] pom;
                pom = mainPopHolder;
                mainPopHolder = trialsPopHolder;
                trialsPopHolder = pom;

                Fmax = this.FindMax(mainPopHolder, fitnessPointers);
                Fmin = this.FindMin(mainPopHolder, fitnessPointers);

                //uCr - computation
                uCr = (1 - c) * uCr + c * meanA(Scr, popEvals, SArchivePointers, N);
                //uF - computation
                uF = (1 - c) * uF + c * meanL(Sf, SArchivePointers, N);

                popEvals += N;
            }

            condition = false;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i].Set();
            }

            outerCountdown.Wait();

            Result result = new Result();

            result.FuncEvals = popEvals;
            result.Fmin = Fmin;
            result.Point = Common.Common.GetBestMember(mainPopHolder, fitnessPointers, N, dim);


            return result;
        }

        private double meanL(double[] Sf, int[] SArchivePointers, int N)
        {
            double f2 = 0;

            for (int i = 0; i < N; i++)
            {
                f2 += (Sf[SArchivePointers[i]] * Sf[SArchivePointers[i]]);
            }

            double f = 0;

            for (int i = 0; i < N; i++)
            {
                f += Sf[SArchivePointers[i]];
            }

            double ret = f2 / f;

            if (double.IsNaN(ret) || double.IsInfinity(ret))
            {
                return 0;
            }
            else
            {
                return ret;
            }
        }

        private double meanA(double[] Scr, int popEvals, int[] SArchivePointers, int N)
        {
            double mean = 0;
            double len = 0;

            for (int i = 0; i < N; i++)
            {
                mean += Scr[SArchivePointers[i]];

                if (Scr[SArchivePointers[i]] != 0)
                {
                    len++;
                }
            }

            if (len == 0 || mean == 0)
            {
                return 0;
            }

            return mean / len;
        }

        private int GetBestPMemberIndex(double[] mainPopHolder, int[] fitnessPointers, double p, int N, RngStream localStream)
        {
            int count = (int)Math.Round(N * p);

            int[] bestPMembersIndicies = new int[count];

            double latestMin = double.MinValue;

            for (int i = 0; i < count; i++)
            {
                double minVal = double.MaxValue;

                for (int j = 0; j < N; j++)
                {
                    if (minVal > mainPopHolder[fitnessPointers[j]] && mainPopHolder[fitnessPointers[j]] > latestMin)
                    {
                        bestPMembersIndicies[i] = j;
                        minVal = mainPopHolder[fitnessPointers[j]];
                    }
                }

                latestMin = minVal;
            }

            return bestPMembersIndicies[localStream.Next(count)];
        }


    }
}
